###Space, Time, Space vs. Time#####
###02/01/2022######################
###Generate Plots##################


#Generic working directory
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

#Get packages, attach them
pks=c("ggplot2", "reshape2")
needed=setdiff(pks,installed.packages())

if( length(needed) > 0){
  options(Ncpus=parallel::detectCores() %/% 1.25)
  install.packages(needed)
}

sapply(pks, library, character.only=T)
session = sessionInfo()


beta.results.static.all <- matrix(NA, nrow=7, ncol=6)
beta.results.ldv.all <- matrix(NA, nrow=7, ncol=6)
beta.results.sar.all <- matrix(NA, nrow=7, ncol=6)
beta.results.stadl.all <- matrix(NA, nrow=7, ncol=6)
phi.results.ldv.all <- matrix(NA, nrow=7, ncol=6)
phi.results.stadl.all <- matrix(NA, nrow=7, ncol=6)
rho.results.sar.all <- matrix(NA, nrow=7, ncol=6)
rho.results.stadl.all <- matrix(NA, nrow=7, ncol=6)
lm.test.all <- matrix(NA, nrow=7, ncol=6) 

load("FINALphi_0.RData")

beta.results.static.all[,1] <- beta.results.static[1,]
beta.results.ldv.all[,1]  <- beta.results.ldv[1,]
beta.results.sar.all[,1]  <- beta.results.sar[1,]
beta.results.stadl.all[,1]  <- beta.results.stadl[1,]
phi.results.ldv.all[,1]  <- phi.results.ldv[1,]
phi.results.stadl.all[,1]  <- phi.results.stadl[1,]
rho.results.sar.all[,1]  <- rho.results.sar[1,]
rho.results.stadl.all[,1]  <- rho.results.stadl[1,]
lm.test.all[,1] <- lm.test[1,]

load("FINALphi_0.1.RData")

beta.results.static.all[,2] <- beta.results.static[1,]
beta.results.ldv.all[,2]  <- beta.results.ldv[1,]
beta.results.sar.all[,2]  <- beta.results.sar[1,]
beta.results.stadl.all[,2]  <- beta.results.stadl[1,]
phi.results.ldv.all[,2]  <- phi.results.ldv[1,]
phi.results.stadl.all[,2]  <- phi.results.stadl[1,]
rho.results.sar.all[,2]  <- rho.results.sar[1,]
rho.results.stadl.all[,2]  <- rho.results.stadl[1,]
lm.test.all[,2] <- lm.test[1,]


load("FINALphi_0.2.RData")

beta.results.static.all[,3] <- beta.results.static[1,]
beta.results.ldv.all[,3]  <- beta.results.ldv[1,]
beta.results.sar.all[,3]  <- beta.results.sar[1,]
beta.results.stadl.all[,3]  <- beta.results.stadl[1,]
phi.results.ldv.all[,3]  <- phi.results.ldv[1,]
phi.results.stadl.all[,3]  <- phi.results.stadl[1,]
rho.results.sar.all[,3]  <- rho.results.sar[1,]
rho.results.stadl.all[,3]  <- rho.results.stadl[1,]
lm.test.all[,3] <- lm.test[1,]


load("FINALphi_0.3.RData")

beta.results.static.all[,4] <- beta.results.static[1,]
beta.results.ldv.all[,4]  <- beta.results.ldv[1,]
beta.results.sar.all[,4]  <- beta.results.sar[1,]
beta.results.stadl.all[,4]  <- beta.results.stadl[1,]
phi.results.ldv.all[,4]  <- phi.results.ldv[1,]
phi.results.stadl.all[,4]  <- phi.results.stadl[1,]
rho.results.sar.all[,4]  <- rho.results.sar[1,]
rho.results.stadl.all[,4]  <- rho.results.stadl[1,]
lm.test.all[,4] <- lm.test[1,]

load("FINALphi_0.4.RData")

beta.results.static.all[,5] <- beta.results.static[1,]
beta.results.ldv.all[,5]  <- beta.results.ldv[1,]
beta.results.sar.all[,5]  <- beta.results.sar[1,]
beta.results.stadl.all[,5]  <- beta.results.stadl[1,]
phi.results.ldv.all[,5]  <- phi.results.ldv[1,]
phi.results.stadl.all[,5]  <- phi.results.stadl[1,]
rho.results.sar.all[,5]  <- rho.results.sar[1,]
rho.results.stadl.all[,5]  <- rho.results.stadl[1,]
lm.test.all[,5] <- lm.test[1,]

load("FINALphi_0.5.RData")

beta.results.static.all[,6] <- beta.results.static[1,]
beta.results.ldv.all[,6]  <- beta.results.ldv[1,]
beta.results.sar.all[,6]  <- beta.results.sar[1,]
beta.results.stadl.all[,6]  <- beta.results.stadl[1,]
phi.results.ldv.all[,6]  <- phi.results.ldv[1,]
phi.results.stadl.all[,6]  <- phi.results.stadl[1,]
rho.results.sar.all[,6]  <- rho.results.sar[1,]
rho.results.stadl.all[,6]  <- rho.results.stadl[1,]
lm.test.all[,6] <- lm.test[1,]


##Figure 6 - LDV phi bias 

rho = seq(0.0, 0.30, 0.05)
phi.results.ldv.all = cbind(phi.results.ldv.all, rho)
colnames(phi.results.ldv.all ) <- c("P-0.00","P-0.10","P-0.20","P-0.30","P-0.40","P-0.50","rho")

cplot_data6 = reshape(as.data.frame(phi.results.ldv.all), 
                     varying = c("P-0.00","P-0.10","P-0.20","P-0.30","P-0.40","P-0.50"), 
                     timevar = "phi",
                     direction = "long",
                     idvar = "rho",  
                     sep="-"
)


cplot_data6$phi <- as.factor(cplot_data6$phi)

p = ggplot(data=cplot_data6, aes(x=rho, y=P, group=phi, linetype=phi)) 
p + geom_line(size = 0.75) +  
  theme_bw(base_size = 15) +
  theme(legend.key.width = unit(1.5,"cm")) +
  ylab(expression(hat(phi)~-phi))+
  xlab( expression(rho[Y])) +
  labs(linetype=expression(phi[Y]))

ggsave("Figure6.pdf", height=4, width=6, units="in", dpi=600)


##Figure 7 - LDV beta bias 
rho = seq(0.0, 0.30, 0.05)
beta.results.ldv.all = cbind(beta.results.ldv.all, rho)
colnames(beta.results.ldv.all ) <- c("P-0.00","P-0.10","P-0.20","P-0.30","P-0.40","P-0.50","rho")

cplot_data7 = reshape(as.data.frame(beta.results.ldv.all), 
                     varying = c("P-0.00","P-0.10","P-0.20","P-0.30","P-0.40","P-0.50"), 
                     timevar = "phi",
                     direction = "long",
                     idvar = "rho",  
                     sep="-"
)


cplot_data7$phi <- as.factor(cplot_data7$phi)

p = ggplot(data=cplot_data7, aes(x=rho, y=P, group=phi, linetype=phi)) 
p + geom_line(size = 0.75) +  
  theme_bw(base_size = 15) +
  theme(legend.key.width = unit(1.5,"cm")) +
  ylab(expression(hat(beta)[LDV]~-beta))+
  xlab( expression(rho[Y])) + 
  labs(linetype=expression(phi[Y]))

ggsave("Figure7.pdf", height=4, width=6, units="in", dpi=600)


##Figure 8 - Response-Path Estimates (direct, indirect, total)


beta.ldv <- beta.results.ldv.all[7,6] + beta 
phi.ldv <- phi.results.ldv.all[7,6] + 0.5

beta.static <- beta.results.static.all[7,6] + beta 

ldv.part= solve(IT - phi.ldv*TL)*beta.ldv #using the average values across the simulations 
true.part = solve(IT - 0.5*TL - 0.3*WT)*2

ldv.effect = matrix(NA, nrow=19, ncol=1)

for(i in 1:19){ 
  ldv.effect[i] = mean(colSums(ldv.part[(1+(50*(i-1))):(50*i),1:50])) 
}

true.direct.effect = matrix(NA, nrow=19, ncol=50)
true.indirect.effect = matrix(NA, nrow=19, ncol=50)
true.tot.effect = matrix(NA, nrow=19, ncol=50)


for(i in 1:19){ 
  true.direct.effect[i,] = diag(true.part[(1+(50*(i-1))):(50*i),1:50])
}

true.direct.mean = rep(rowMeans(true.direct.effect),50)

for(i in 1:19){ 
  true.indirect.effect[i,] = colSums(true.part[(1+(50*(i-1))):(50*i),1:50]) - diag(true.part[(1+(50*(i-1))):(50*i),1:50])
}

true.indirect.mean = rep(rowMeans(true.indirect.effect),50)


for(i in 1:19){ 
  true.tot.effect[i,] = colSums(true.part[(1+(50*(i-1))):(50*i),1:50])
}

true.tot.mean = rep(rowMeans(true.tot.effect),50)


year = seq(1, 19, 1)
lrrp.direct.effects = cbind(true.direct.effect,year)
lrrp.indirect.effects = cbind(true.indirect.effect,year)
lrrp.tot.effects = cbind(true.tot.effect,year)

lrrp.direct.effects.long = melt(lrrp.direct.effects, 'year')
lrrp.direct.effects.long = lrrp.direct.effects.long[1:950,]

lrrp.indirect.effects.long = melt(lrrp.indirect.effects, 'year')
lrrp.indirect.effects.long = lrrp.indirect.effects.long[1:950,]

lrrp.tot.effects.long = melt(lrrp.tot.effects, 'year')
lrrp.tot.effects.long = lrrp.tot.effects.long[1:950,]



unit = rep(1:50, each=19)
ldv.effect = rep(ldv.effect,50)


static.effect.t = c(beta.static,rep(0,18))
static.effect = rep(static.effect.t,50)

data_fig8a = cbind(lrrp.direct.effects.long,ldv.effect,true.direct.mean,static.effect,unit)
data_fig8a$unit = as.factor(data_fig8a$unit)

ggplot(data=data_fig8a, aes(x=year,y=value,group=unit)) + 
  theme_bw()+
  ylab(expression('Response of y'))+
  xlab("Period") +
  geom_line(aes(x=year, y=value, group=unit), colour="grey") +
  geom_line(aes(x=year, y=ldv.effect), colour="black", linetype = "dashed",size=1.0) +   
  geom_line(aes(x=year, y=true.direct.mean), colour="black",size=0.75) + 
  geom_line(aes(x=year, y=static.effect), colour="black", linetype="dotted", size=0.75) + 
  scale_x_continuous(breaks=c(1,6,11,16),
                   labels=c("t", "t+5","t+10","t+15"))


ggsave("Figure8a.pdf", height=4, width=6, units="in", dpi=600)


ldv.indirect.effect <- rep(0,950)
static.indirect.effect <- rep(0,950)

data_fig8b = cbind(lrrp.indirect.effects.long,ldv.indirect.effect,static.indirect.effect,true.indirect.mean,unit)
data_fig8b$unit = as.factor(data_fig8b$unit)

ggplot(data=data_fig8b, aes(x=year,y=value,group=unit)) + 
  theme_bw()+
  ylab(expression('Response of y'))+
  xlab("Period") +
  geom_line(aes(x=year, y=value, group=unit), colour="grey") +
  geom_line(aes(x=year, y=ldv.indirect.effect), colour="black", linetype = "dashed",size=1.0) +   
  geom_line(aes(x=year, y=true.indirect.mean), colour="black",size=0.75) + 
  geom_line(aes(x=year, y=static.indirect.effect), colour="black", linetype="dotted", size=0.75) + 
  scale_x_continuous(breaks=c(1,6,11,16),
                     labels=c("t", "t+5","t+10","t+15"))

ggsave("Figure8b.pdf", height=4, width=6, units="in", dpi=600)



data_fig8c = cbind(lrrp.tot.effects.long,ldv.effect,static.effect,true.tot.mean,unit)
data_fig8c$unit = as.factor(data_fig8c$unit)

ggplot(data=data_fig8c, aes(x=year,y=value,group=unit)) + 
  theme_bw()+
  ylab(expression('Response of y'))+
  xlab("Period") +
  geom_line(aes(x=year, y=value, group=unit), colour="grey") +
  geom_line(aes(x=year, y=ldv.effect), colour="black", linetype = "dashed",size=1.0) +   
  geom_line(aes(x=year, y=true.tot.mean), colour="black",size=0.75) +
  geom_line(aes(x=year, y=static.effect), colour="black", linetype="dotted", size=0.75) + 
  scale_x_continuous(breaks=c(1,6,11,16),
                     labels=c("t", "t+5","t+10","t+15"))

ggsave("Figure8c.pdf", height=4, width=6, units="in", dpi=600)


##Figure 9 - Cumulative Response Path 

cum.ldv.effect =  cumsum(ldv.effect[1:19])
cum.static.effect = cumsum(static.effect[1:19])
cum.direct.effect = apply(true.direct.effect, 2, cumsum)
cum.indirect.effect = apply(true.indirect.effect, 2, cumsum)
cum.tot.effect = apply(true.tot.effect, 2, cumsum)



cum.direct.mean = rowMeans(cum.direct.effect)
cum.indirect.mean = rowMeans(cum.indirect.effect)
cum.tot.mean = rowMeans(cum.tot.effect)

cum.direct.effect2 = cbind(cum.direct.effect, year)
cum.indirect.effect2 = cbind(cum.indirect.effect, year)
cum.tot.effect2 = cbind(cum.tot.effect, year)


cum.direct.effect.long = melt(cum.direct.effect2, 'year')
cum.direct.effect.long = cum.direct.effect.long[1:950,]

cum.indirect.effect.long = melt(cum.indirect.effect2, 'year')
cum.indirect.effect.long = cum.indirect.effect.long[1:950,]

cum.tot.effect.long = melt(cum.tot.effect2, 'year')
cum.tot.effect.long = cum.tot.effect.long[1:950,]


unit = rep(1:50, each=19)

cum.ldv.effect.long = rep(cum.ldv.effect,50)
cum.static.effect.long = rep(cum.static.effect,50)
cum.direct.mean.long = rep(cum.direct.mean,50)
cum.indirect.mean.long = rep(cum.indirect.mean,50)
cum.tot.mean.long = rep(cum.tot.mean,50)


data_fig9a = cbind(cum.direct.effect.long,cum.ldv.effect.long,cum.static.effect.long,cum.direct.mean,unit)
data_fig9a$unit = as.factor(data_fig9a$unit)

ggplot(data=data_fig9a, aes(x=year,y=value,group=unit)) + 
  theme_bw()+
  ylab(expression('Response of y'))+
  xlab("Period") +
  geom_line(aes(x=year, y=value, group=unit), colour="grey") +
  geom_line(aes(x=year, y=cum.ldv.effect.long), colour="black", linetype = "dashed",size=1.0) +   
  geom_line(aes(x=year, y=cum.direct.mean), colour="black",size=0.75) + 
  geom_line(aes(x=year, y=cum.static.effect.long), colour="black", linetype="dotted", size=0.75) + 
  scale_x_continuous(breaks=c(1,6,11,16),
                     labels=c("t", "t+5","t+10","t+15"))

ggsave("Figure9a.pdf", height=4, width=6, units="in", dpi=600)

cum.ldv.indirect.effect <- rep(0,950)
cum.static.indirect.effect <- rep(0,950)

data_fig9b = cbind(cum.indirect.effect.long,cum.ldv.indirect.effect,cum.static.indirect.effect,cum.indirect.mean,unit)
data_fig9b$unit = as.factor(data_fig9b$unit)

ggplot(data=data_fig9b, aes(x=year,y=value,group=unit)) + 
  theme_bw()+
  ylab(expression('Response of y'))+
  xlab("Period") +
  geom_line(aes(x=year, y=value, group=unit), colour="grey") +
  geom_line(aes(x=year, y=cum.ldv.indirect.effect), colour="black", linetype = "dashed",size=1.0) +   
  geom_line(aes(x=year, y=cum.indirect.mean), colour="black",size=0.75) +
  geom_line(aes(x=year, y=cum.static.indirect.effect), colour="black", linetype="dotted", size=0.75) + 
  scale_x_continuous(breaks=c(1,6,11,16),
                     labels=c("t", "t+5","t+10","t+15"))

ggsave("Figure9b.pdf", height=4, width=6, units="in", dpi=600)


data_fig9c = cbind(cum.tot.effect.long,cum.ldv.effect.long,cum.static.effect.long,cum.tot.mean,unit)
data_fig9c$unit = as.factor(data_fig9c$unit)

ggplot(data=data_fig9c, aes(x=year,y=value,group=unit)) + 
  theme_bw()+
  ylab(expression('Response of y'))+
  xlab("Period") +
  geom_line(aes(x=year, y=value, group=unit), colour="grey") +
  geom_line(aes(x=year, y=cum.ldv.effect.long), colour="black", linetype = "dashed",size=1.0) +   
  geom_line(aes(x=year, y=cum.tot.mean), colour="black",size=0.75) +
  geom_line(aes(x=year, y=cum.static.effect.long), colour="black", linetype="dotted", size=0.75) + 
  scale_x_continuous(breaks=c(1,6,11,16),
                     labels=c("t", "t+5","t+10","t+15"))

ggsave("Figure9c.pdf", height=4, width=6, units="in", dpi=600)



##Figure 10 - SAR rho bias 


rho = seq(0.0, 0.30, 0.05)
rho.results.sar.all = cbind(rho.results.sar.all, rho)
colnames(rho.results.sar.all ) <- c("P-0.00","P-0.10","P-0.20","P-0.30","P-0.40","P-0.50","rho")

cplot_data10 = reshape(as.data.frame(rho.results.sar.all), 
                      varying = c("P-0.00","P-0.10","P-0.20","P-0.30","P-0.40","P-0.50"), 
                      timevar = "phi",
                      direction = "long",
                      idvar = "rho",  
                      sep="-"
)


cplot_data10$phi <- as.factor(cplot_data10$phi)

p = ggplot(data=cplot_data10, aes(x=rho, y=P, group=phi, linetype=phi)) 
p + geom_line(size = 0.75) +  
  theme_bw(base_size = 15) +
  theme(legend.key.width = unit(1.5,"cm")) +
  ylab(expression(hat(rho)~-rho))+
  xlab( expression(rho[Y])) + 
  labs(linetype=expression(phi[Y]))

ggsave("Figure10.pdf", height=4, width=6, units="in", dpi=600)


##Figure 11 - SAR beta bias 

rho = seq(0.0, 0.30, 0.05)
beta.results.sar.all = cbind(beta.results.sar.all, rho)
colnames(beta.results.sar.all ) <- c("P-0.00","P-0.10","P-0.20","P-0.30","P-0.40","P-0.50","rho")

cplot_data11 = reshape(as.data.frame(beta.results.sar.all), 
                      varying = c("P-0.00","P-0.10","P-0.20","P-0.30","P-0.40","P-0.50"), 
                      timevar = "phi",
                      direction = "long",
                      idvar = "rho",  
                      sep="-"
)


cplot_data11$phi <- as.factor(cplot_data11$phi)

p = ggplot(data=cplot_data11, aes(x=rho, y=P, group=phi, linetype=phi)) 
p + geom_line(size = 0.75) +  
  theme_bw(base_size = 15) +
  theme(legend.key.width = unit(1.5,"cm")) +
  ylab(expression(hat(beta)[SAR]~-beta))+
  xlab( expression(rho[Y])) + 
  labs(linetype=expression(phi[Y]))


ggsave("Figure11.pdf", height=4, width=6, units="in", dpi=600)




##Figure A1 - Power LM test 


rho = seq(0.0, 0.30, 0.05)
lm.test.all = cbind(lm.test.all, rho)
colnames(lm.test.all ) <- c("P-0.00","P-0.10","P-0.20","P-0.30","P-0.40","P-0.50","rho")

cplot_dataA1 = reshape(as.data.frame(lm.test.all), 
                      varying = c("P-0.00","P-0.10","P-0.20","P-0.30","P-0.40","P-0.50"), 
                      timevar = "phi",
                      direction = "long",
                      idvar = "rho",  
                      sep="-"
)


cplot_dataA1$phi <- as.factor(cplot_dataA1$phi)

p = ggplot(data=cplot_dataA1, aes(x=rho, y=P/tr, group=phi, linetype=phi)) 
p + geom_line() +  
  theme_bw(base_size = 20) +
  ylab(expression(False~positive~rate))+
  xlab( expression(rho[Y])) + 
  labs(linetype=expression(phi[Y]))



ggsave("FigureA1.pdf", height=4, width=6, units="in", dpi=600)


##Figure A2 - STADL Beta 

rho = seq(0.0, 0.30, 0.05)
beta.results.stadl.all = cbind(beta.results.stadl.all, rho)
colnames(beta.results.stadl.all ) <- c("P-0.00","P-0.10","P-0.20","P-0.30","P-0.40","P-0.50","rho")

cplot_dataA2 = reshape(as.data.frame(beta.results.stadl.all), 
                      varying = c("P-0.00","P-0.10","P-0.20","P-0.30","P-0.40","P-0.50"), 
                      timevar = "phi",
                      direction = "long",
                      idvar = "rho",  
                      sep="-"
)


cplot_dataA2$phi <- as.factor(cplot_dataA2$phi)

p = ggplot(data=cplot_dataA2, aes(x=rho, y=P, group=phi, linetype=phi)) 
p + geom_line() +  
  ylim(0.1,-0.1) +
  theme_bw(base_size = 20) +
  ylab(expression(hat(beta)[STADL]~-beta))+
  xlab( expression(rho[Y])) + 
  labs(linetype=expression(phi[Y]))

ggsave("FigureA2.pdf", height=4, width=6, units="in", dpi=600)


##Figure A3 - STADL Rho 

rho = seq(0.0, 0.30, 0.05)
rho.results.stadl.all = cbind(rho.results.stadl.all, rho)
colnames(rho.results.stadl.all ) <- c("P-0.00","P-0.10","P-0.20","P-0.30","P-0.40","P-0.50","rho")

cplot_dataA3 = reshape(as.data.frame(rho.results.stadl.all), 
                      varying = c("P-0.00","P-0.10","P-0.20","P-0.30","P-0.40","P-0.50"), 
                      timevar = "phi",
                      direction = "long",
                      idvar = "rho",  
                      sep="-"
)


cplot_dataA3$phi <- as.factor(cplot_dataA3$phi)

p = ggplot(data=cplot_dataA3, aes(x=rho, y=P, group=phi, linetype=phi)) 
p + geom_line() +  
  theme_bw(base_size = 20) +
  ylim(0.1,-0.1) +
  ylab(expression(hat(rho)[STADL]~-rho))+
  xlab( expression(rho[Y])) + 
  labs(linetype=expression(phi[Y]))

ggsave("FigureA3.pdf", height=4, width=6, units="in", dpi=600)

##Figure A4 - STADL Phi 


rho = seq(0.0, 0.30, 0.05)
phi.results.stadl.all = cbind(phi.results.stadl.all, rho)
colnames(phi.results.stadl.all ) <- c("P-0.00","P-0.10","P-0.20","P-0.30","P-0.40","P-0.50","rho")

cplot_dataA4 = reshape(as.data.frame(phi.results.stadl.all), 
                      varying = c("P-0.00","P-0.10","P-0.20","P-0.30","P-0.40","P-0.50"), 
                      timevar = "phi",
                      direction = "long",
                      idvar = "rho",  
                      sep="-"
)


cplot_dataA4$phi <- as.factor(cplot_dataA4$phi)

p = ggplot(data=cplot_dataA4, aes(x=rho, y=P, group=phi, linetype=phi)) 
p + geom_line() +  
  theme_bw(base_size = 20) +
  ylim(0.1,-0.1) +
  ylab(expression(hat(phi)[STADL]~-phi))+
  xlab( expression(rho[Y])) + 
  labs(linetype=expression(phi[Y]))

ggsave("FigureA4.pdf", height=4, width=6, units="in", dpi=600)


##Figure A5 - Indirect Response Path 

#beta.temp <- beta.results.ldv.all[16,6] + beta 
#phi.temp <- phi.results.ldv.all[16,6] + 0.5


true.part = solve(IT - 0.5*TL - 0.3*WT)*2
true.effect = matrix(NA, nrow=19, ncol=50)


for(i in 1:19){ 
  true.effect[i,] = colSums(true.part[(1+(50*(i-1))):(50*i),1:50]) - diag(true.part[(1+(50*(i-1))):(50*i),1:50])
}

true.mean = rep(rowMeans(true.effect),50)

year = seq(1, 19, 1)
lrrp.effects = cbind(true.effect, year)


data_figA5 = melt(lrrp.effects, 'year')
data_figA5 = data_figA5[1:950,]
unit = rep(1:50, each=19)

data_figA5 = cbind(data_figA5,true.mean, unit)
data_figA5$unit = as.factor(data_figA5$unit)

ggplot(data=data_figA5, aes(x=year,y=value,group=unit)) + 
  theme_bw()+
  ylab(expression('Indirect effect of change in X'[i]*' in year t'))+
  xlab("Years after t") +
  geom_line(aes(x=year, y=value, group=unit), colour="grey") +
  geom_line(aes(x=year, y=true.mean), colour="black",size=1.0) 


ggsave("FigureA5.pdf", height=4, width=6, units="in", dpi=600)


